load nursery.mat;
reset(RandStream.getDefaultStream)
nurseryTest = nursery(randperm(size(nursery,1)),:);
train = nurseryTest(1:size(nursery,1)/2,:);
test = nurseryTest(size(nursery,1)/2+1:end,:);

% P(x1,x2,...,xn,Y) = P(x1=1|Y=1)P(x2=5|Y=1)...P(xn=4|Y=1)P(Y=1)
%P(x1|Y)=  multinomial/multinoulli distribution

numXValues = size(nurseryTest,2) -1;

%find P(Y)
yFreq = zeros(5,1);

for trainIndex=1:5
    yFreq(trainIndex,1) = sum(train(:,numXValues+1)==trainIndex);
end
yTheta = (yFreq+ones(5,1))./(size(train,1)+max(train(:,numXValues+1)));

%trainResults has dimensions d BY possible values of x BY possible values of Y
trainResults = zeros(numXValues,5,5);

for xIndex=1:numXValues
    for xvalue=1:5
        for yvalue=1:5
            trainResults(xIndex,xvalue,yvalue) = (sum(train(:,numXValues+1)==yvalue&train(:,xIndex)==xvalue)+1) / (sum(train(:,numXValues+1)==yvalue)+max(train(:,xIndex)));
        end
    end
end

testProbabilities =  ones(size(test,1),5);
for testIndex = 1:size(test,1)
    for yIndex = 1:5
        for xIndex = 1:numXValues        
            testProbabilities(testIndex,yIndex) = testProbabilities(testIndex,yIndex) * trainResults(xIndex,test(testIndex,xIndex),yIndex);
        end
        testProbabilities(testIndex,yIndex) = testProbabilities(testIndex,yIndex) * yTheta(yIndex,1);
    end
end

[~,testIndexResults] = max(testProbabilities');
testIndexResults = testIndexResults';

A = [testIndexResults test(:,numXValues+1)];
accuracy = sum(A(:,1)==A(:,2)) / size(A,1);

logLikelihood = 0;
for index=1:size(train,1)
    for dimensionIndex=1:numXValues
        logLikelihood = logLikelihood + log(trainResults(dimensionIndex,train(index, dimensionIndex),train(index, numXValues + 1)));
    end
    logLikelihood = logLikelihood + log(yTheta(train(index, numXValues + 1)));
end

testLogLikelihood = 0;
for index=1:size(test,1)
    for dimensionIndex=1:numXValues
        testLogLikelihood = testLogLikelihood + log(trainResults(dimensionIndex,test(index, dimensionIndex),test(index, numXValues + 1)));
    end
    testLogLikelihood = testLogLikelihood + log(yTheta(test(index, numXValues + 1)));
end
